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Introduction 


Conventional  mammography  is  limited  in  its  sensitivity  for  detecting  subtle  breast 
pathological  changes,  since  the  imaging  relies  on  the  small  differences  in  x-ray  attenuation 
between  the  lesions  and  breast  tissues.  As  x-ray  traverses  a  breast,  it  gets  not  only  attenuated, 
but  also  phase-shifted.  Utilizing  the  highly  sensitive  x-ray  phase-shifts,  the  phase  sensitive  x- 
ray  imaging  has  the  potential  of  greatly  increasing  x-ray  imaging  sensitivity  and  specificity  and 
reducing  radiation  doses  associated  with  the  imaging. 

The  long  term  objective  of  the  project  is  to  develop  a  low-dose  and  quantitative  x-ray 
phase  imaging  technique  for  facilitating  breast  cancer  detection.  The  Specific  Aims  of  the 
project  of  the  three-year  period  are:  (1)  Develop  a  prototype  phase-imaging  system  enabling 
the  phase  retrieval,  that  is,  the  reconstruction  of  objects  phase-maps.  The  system  hardware 
comprises  a  micro-focus  tube  operating  at  high  tube  voltages,  a  high  resolution  photostimulable 
phosphor  plate  (CR-plate)  based  detector  system.  The  core  algorithms  for  breast  phase-map 
reconstruction  will  be  developed  to  retrieve  a  breast  phase  map  from  a  single  recorded  image. 
(2)  Validate  the  accuracy  of  the  reconstructed  tissue  projected  electron  densities;  validate  the 
many-fold  radiation  dose  reduction  achieved  with  the  proposed  system;  conduct  subjective 
measurements  to  characterize  the  performance  of  the  proposed  system. 


Body 


In  this  period  we  have  performed  following  studies:  (A)  Since  the  robustness  of  phase 
retrieval  methods  is  of  critical  importance  for  limiting  and  reducing  radiation  doses  involved  in  x- 
ray  phase  imaging,  in  this  period  we  have  compared  the  robustness  of  the  attenuation-partition 
based  (AP-based)  phase  retrieval  method  that  we  developed  to  robustness  of  the  Transport  of 
Intensity  Equation  based  (TIE-based)  method.  (B)  In  order  to  develop  low-dose  and  quantitative 
x-ray  phase  imaging,  in  this  period  we  have  conducted  phantom  imaging  studies  with  100  kVp- 
140  kVp  tube  voltages  and  performed  the  phase  retrievals  to  retrieve  phantom  phase  maps  for 
fully  exhibiting  the  phase  contrast  and  providing  quantitative  information  on  the  phantom 
electron  densities. 


(A).  Experimental  performance  evaluation  of  the  attenuation-partition  based  iterative 
phase  retrieval  method 


Introduction:  In  the  Period  2  we  developed  a  novel  phase  retrieval  method,  the 
attenuation-partition  based  (AP-based)  iterative  phase  retrieval  method.  The  phase  retrieval  is  a 
procedure  of  retrieving  the  phase-shift  map^3(r)  of  an  object  from  its  phase-sensitive 

projections,  where  ^(F)  = -;irg  J  p^{r ,s)ds  yOg(F, 5)  denotes  the  tissue  electron  density,  = 

ray 

2.81 8  X  1 0'^®  m  is  the  classical  electron  radius  and  X  is  x-ray  wavelength,  the  integral  is 
computed  over  the  ray  path.  A  retrieved  phase  map  is  able  to  fully  exhibit  tissue  phase  contrast 
and  provide  a  quantitative  map  of  the  object's  projected  electron  densities,  which  could  be  used 
for  quantitative  tissue  characterizations.  Moreover,  performing  phase  retrieval  is  necessary  for 
reconstructing  volumetric  3-D  maps  of  tissue  attenuation  coefficients  and  refractive  indices 
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respectively,  and  for  eliminating  the  phase-contrast  caused  artifacts  in  the  volumetric  3-D 
images. 

Examining  the  robustness  of  phase  retrieval  methods  against  the  projection  noise,  such 
as  the  x-ray  photon  quantum  noise  involved  in  the  projection  acquisitions,  is  an  important  task  in 
developing  phase  imaging  techniques  for  future  medical  imaging  applications.  If  a  phase 
retrieval  method  is  not  robust  against  the  noise,  the  method  could  become  unstable  and  fail  for 
phase  retrieval  if  the  acquired  projection  images  were  contaminated  with  substantial  noise.  On 
the  other  hand,  suppressing  x-ray  quantum  noise  may  require  using  high  radiation  doses  in  the 
image  acquisitions.  For  clinical  applications  it  is  imperative  to  limit  and  reduce  radiation  doses 
involved,  hence  it  is  critical  to  develop  robust  phase  retrieval  methods  for  future  clinical 
applications  of  phase  contrast  imaging.  The  most  commonly  used  method  for  phase  retrieval  in 
the  literature  is  based  on  the  Transport  of  Intensity  Equation  (TIE)  for  phase  contrast 
projections  [Allen  et  al.  2001].  While  the  TIE-based  phase  retrieval  method  is  computational 
effective  for  the  cases  with  strong  phase  contrast  effects  and  low  noise  levels,  we  found  that  it 
would  fail  for  phase  retrieval  when  the  acquired  projection  images  were  contaminated  with 
relatively  strong  noises  as  in  clinical  applications  where  radiation  dose  constraints  dictates 
relatively  high  noise  levels.  To  improve  the  robustness  of  phase  retrieval,  we  recently 
developed  a  novel  method  called  the  attenuation  partition  (AP)  based  method,  or  the  AP-based 
method  for  short.  In  this  period  we  compared  the  robustness  of  two  phase-retrieval  methods,  as 
the  TIE-based  method  vs.  the  AP-based  method,  by  analyzing  the  phase  maps  retrieved  from 
the  experimental  images  of  a  phantom. 


Method:  We  applied  these  two  methods  to  experimental  projection  images  of  an  air- 
bubble  wrap  phantom  for  retrieving  the  phase  map  of  the  bubble  wrap.  According  to  the  TIE- 
based  method,  one  acquires  two  images:  one  contact  radiograph,  and  one  phase  contrast 
image.  The  phase  map  (p{r)  is  retrieved  as 


^ (r )  =  - [inMlXR^ )  j V  ■  [ V ( 


[Allen  et  al.  2001],  where  aI  denotes 


the  its  x-ray  attenuation  map  of  an  object  obtained  from  its  contact  radiograph,  A  is  the  average 
x-ray  wavelength,  I  the  phase-sensitive  projection  intensity,  /,„the  entrance  x-ray  intensity, 
the  source-object  distance  set  in  the  projection,  is  the  object-detector  distance  and 
M  =  (7?j  +7?2 )/./?!  the  magnification  factor  in  the  projection.  In  above  formula  the  operator  V 


denotes  the  two-dimensional  transverse  gradient  differential  operator,  and  V  ^denotes  the 
inverse  Laplacian  operator.  We  noted  that  the  Achilles  heel  of  the  TIE-based  method  lies  at  the 

inverse  Laplacian  operator  involved  in  the  method.  This  singularity  may  amplify  the  noise 
randomly  embedded  in  the  acquired  images  and  result  in  instability  in  phase  retrievals.  To  mend 
the  instability  of  the  TIE-method,  one  may  use  the  Tikhonov  regularization  technique,  which 

regularizes  the  inverse  Laplacian  operator  attempting  to  tame  the  singularity-caused 
problems,  but  its  effectiveness  is  very  limited  [Allen  et  al.  2001], 

In  order  to  get  rid  of  the  singular  pseudo-differential  operator  involved  in  the  TIE- 
based  phase  retrievals,  we  recently  developed  a  novel  phase  retrieval  method:  the  AP  method. 
Our  idea  is  to  utilize  the  correlation  between  the  x-ray  phase  shift  and  its  attenuation  to 
eliminate  the  singularities  involved  in  phase  retrievals.  According  to  this  the  AP-based  method, 
one  first  applies  the  duality  transform  2)  to  the  acquired  phase  contrast  image /(r)  to  obtain  an 

estimate  of  both  the  Compton  scattering-generated  attenuation  ^^^and  phase  map^(r)of  the 
sample.  Here  the  duality  transform  D  is  implemented  by  using  the  pseudo-differential  operator 
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where 

denotes  the  average  Compton  scattering  cross-section,  the 
classical  electron  radius  the  phase-attenuation  duality  [Wu  et  al. 
2005,  ].  One  then  compares  against  the  acquired 

normalized  radiograph  Al  and  computes  their  differences  as  the 
errors  A  .  Applying  the  Fresnel  diffraction  transform  ^^^.^to  the 
fictitious  transmission  function  exp ,  one  computes  a 
better  estimates  of/^^,  the  projection  intensity  determined  by 

using  the  estimated  object's  electron  densities  and  counting  only 
the  Compton  scattering-generated  attenuation.  With  this  new 
estimate  of  /^^^one  can  start  a  new  round  of  the  iteration  to 

compute  a  new  improved  estimates  of  both  the^^^  and  phase 
map^(r)  of  the  sample.  The  iteration  converges  when  the 

does  not  change  substantially  with  further  iteration  steps. 

In  order  to  compare  the  robustness  of  the  AP-based 
method  against  the  TIE-based  method,  we  applied  the  two 
methods  for  retrieving  the  phase-map  of  a  piece  of  air  bubble 
wrap.  The  x-ray  tube  employed  was  with  a  7- //m  focal  spot.  A  contact  radiograph  and  a  phase 
contrast  image  were  all  acquired  with  40  kVp,  6  mAs  and  a  source-detector  distance  of  1.78  m. 
However,  the  phase  contrast  image  was  acquired  with  a  large  sample-detector  distance  of  1.15 
m  for  allowing  the  exiting  x-rays  to  freely  diffract  over  a  long  distance  and  form  the  interference 
fringes.  The  imaging  detector  employed  was  an  aSe-based  flat-panel  detector  with  a  pixel  pitch 
of  140  jum .  Using  the  acquired  radiograph  and  phase  contrast  image,  the  phase  maps  of  the 

bubble  wrap  were  retrieved  with  the  two  methods,  as  shown  in 
Fig.  1  (with  the  TIE  based  method)  and  Fig.  2  (with  the  AP- 
based  method.),  respectively.  The  image  qualities  and  the 
phase-values  profiles  of  the  two  retrieved  phase  maps  were 
analyzed  for  examining  the  robustness  of  the  two  methods. 

Result:  Fig.  1 -(upper)  shows  the  retrieved  sample's 
phase  map  using  TIE-based  method  .  Apparently  the  retrieval 
with  the  TIE-based  method  is  very  unsatisfactory,  as  the  bubble 
rims  are  hardly  visible  in  the  extremely  cluttered  background  in 
the  phase  map.  In  a  stark  contrast  to  the  TIE-based  method, 
the  bubble  rims  are  prominently  depicted  in  Fig.  2-(upper),  which 
was  retrieved  using  the  AP-method.  In  order  to  gauge  the 
accuracies  of  the  retrieved  phase  values,  we  measured  the 
thicknesses  of  the  flat  bases  of  the  wrap  using  a  caliber  ruler, 
and  we  found  that  the  base  thickness  was  about  0.055  mm. 

The  wrap  base  is  made  of  the  low-density  polyethylene  films. 
According  to  the  film's  molecular  formula  and  its  mass 

density  is  0.925  g/cm^ ,  we  found  that  the  approximate  phase- 
shifts  at  the  flat  bases  of  the  wrap  is  about  -4.2  radian.  On  the 
other  hand,  according  to  the  phase  profile  in  Fig.  2(lower)  with 


Fig.  2.  (Upper)  Retrieved  phase 
maps  of  the  bubble  wrap  with 
the  AP-based  method.  (Lower) 
Profile  of  the  retrieved  phase 
values  along  the  marked  line. 


Fig.  1.  (Upper)  Retrieved  phase 
maps  of  the  bubble  wrap  with 
the  TIE-based  method  and 
mended  with  the  Tikhonov 
regularization.  (Lower)  Profile 
of  the  retrieved  phase  values 
along  the  marked  line. 
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the  AP-based  method,  the  average  phase  shift  at  the  flat  bases  of  the  wrap  is  -1 .05  radian,  the 
two  phase-shift  estimates  are  reasonably  close.  In  comparison,  the  phase  profile  in  Fig. 

1  (lower)  with  the  TIE-based  method,  depicts  just  messy  up-down  peaks  buried  in  cluttered  and 
noisy  background.  According  to  this  profile,  the  phase-shift  values  at  the  bases  fluctuate  over  a 
wide  range  from  0  to  +87.5  radian  with  an  average  of  +27.3  radian,  contradicting  to  the  fact  that 
x-ray  phase  shifts  should  be  negative. 

In  summary,  we  showed  that  the  TIE-based  method  failed  in  retrieving  the  air  bubble 
wrap’s  phase  map.  In  contrast,  the  AP-based  method  retrieved  the  wrap  phase  map 
successfully.  The  retrieved  phase  values  with  the  AP-based  are  reasonably  close  to  the 
estimate  based  on  the  wrap  thickness  measurement.  The  stark  performance  differences  of  the 
two  methods  are  rooted  in  their  different  techniques  employed  to  deal  with  the  singularity 
problem  in  phase  retrievals.  This  comparison  shows  that  the  conventional  TIE-based  phase 
retrieval  method  is  unstable  against  the  noise  in  the  wrap’s  projection  images,  while  the  AP- 
based  phase  retrieval  method  is  shown  in  the  experiment  to  be  superior  to  the  TIE-based 
method  for  the  robustness  in  performing  the  phase  retrieval. 

In  addition,  extending  the  research  work  on  the  phase  retrieval  methods,  as  a  by-product 
of  this  project,  we  made  a  progress  in  understanding  the  apparent  linear  attenuation  coefficients 
in  the  inline  phase  contrast  tomography.  In  the  inline  phase  contrast  x-ray  tomography  the 
reconstructed  apparent  linear  attenuation  coefficient  values  in  the  tomography  may  be  greatly 
larger  than  sample's  linear  attenuation  coefficients  or  even  be  negative.  These  "artifacts"  may 
cause  faulty  interpretation  of  sample  structures,  and  impede  even  qualitative  characterizations 
for  tissues  and  materials.  In  this  period  we  derived  a  general  formula  to  quantitatively  relate  the 
apparent  linear  attenuation  coefficient  values  in  cone-beam  phase  contrast  tomography  to 
sample's  linear  attenuation  coefficients  and  refractive  indices  [Yan  and  Wu  201 1].  This  formula 
overcomes  the  gross  inaccuracy  of  the  existing  formula  in  the  literature  in  analyzing  high- 
resolution  phase  contrast  tomography,  and  it  will  be  useful  for  correctly  interpreting  and 
quantifying  the  apparent  linear  attenuation  coefficients  in  cone-beam  x-ray  phase  contrast 
tomography  in  the  biomedical  and  material  science  applications. 


(B).  Develop  low-dose  and  quantitative  x-ray  phase  imaging  techniques 


Fig.  3.  Phase  contrast  image  of  a 
contrast-detail  phantom  acquired  with 
120  kVp  andM=2. 


Introduction;  As  conventional  mammography 
(digital  or  screen/film)  employs  low  tube  voltages  (25-32 
kVp)  for  breast  imaging,  we  conducted  phantom  imaging 
studies  with  100  kVp-140  kVp  tube  voltages  in  this  period 
for  developing  low-dose  and  quantitative  x-ray  phase 
imaging  techniques.  Our  strategy  is  to  utilize  the  phase 
contrast  to  compensate  for  the  attenuation-contrast  loss 
associated  with  the  high  tube-voltages,  and  to  accomplish 
the  radiation  doses  reduction  with  the  high-voltage  phase 
imaging. 

Method:  The  prototype  system  hardware  consists 
of  a  micro-focus  x-ray  tube  (Hamamatsu  L8121-03)  with  a 
focal  spot  of  50  pm  for  75  W  power  loading,  and  operating 
at  voltages  up  to  150  kVp,  and  a  high-resolution  (43.75 
urn)  CR  system  (Konica  Minolta).  We  have  conducted 
preliminary  imaging  experiments  using  100,  120,  140  kVp 
and  with  four  different  phantoms,  the  ACR  MAP  phantom. 
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a  contrast-detail  phantom,  an  acrylic  step-edge,  and  a  breast  tissue-equivalent  phantom.  We 
images  the  phantoms  with  different  geometric  configurations  including  the  settings  with 
magnification  factors  of  1, 2,  2.5,  3.  Observer  studies  were  performed  utilizing  the  ACR  MAP 
and  contrast-detail  phantoms  for  quantitative  comparisons.  In  this  way  we  seek  to  validate  the 
optimal  geometric  configurations  for  phase  contrast  visibility  as  determined  in  the  first  period 
based  on  the  theoretical  phase  contrast  visibility  analysis  [Wu  et  al.,  2004]. 

In  order  to  accurately  estimate  the  average  glandular  doses  involved  in  the  experiments, 
we  extended  our  previous  tables  of  the  normalized  average  glandular  doses  for  the 

conventional  mammography  to  the values  for  the  high  energy  x-rays,  by  performing  the 

Monte  Carlo  simulations  of  x-ray  photon  transport  and  energy  deposition  in  breast  tissues. 

In  order  to  fully  exhibit  tissue  phase  contrast  and  provide  a  quantitative  information  of  the 
imaged  objects,  performing  phase  retrieval  is  necessary.  While  the  AP-based  phase  retrieval 
method  discussed  in  section  (A)  is  a  general  method,  it  needs  at  least  two  projection  images  for 
phase  retrievals.  In  order  to  further  reduce  radiation  doses  involved,  a  phase  retrieval  method 
requires  only  one  single  phase-sensitive  projection  would  be  even  more  desirable.  Considering 
the  facts  that  breast  consists  of  dominantly  soft  tissues  and  the  x-rays  from  tubes  operating  at 
120-150  kVp  are  dominantly  high-energy  x-rays  (60  keV  or  higher),  we  employed  a  novel 
method  to  retrieve  the  phase  map  of  an  object  by  using  just  one  single  phase  contrast  projection 
image  jWu  et  al.  2005,  Wu  et  al.  2009j.  According  to  this  method,  the  phase  map  of  the  object 
can  be  accurately  retrieved  by  using  the  duality  transform  as  ^(r)  =  (;irg/c7j^^)ln®(/(r^)), 

where  is  the  duality  transform,  as  is  defined  earlier  in  Section  (A),  of  a  phase  contrast 

projection /(?£,)  of  the  object 

Results:  Shown  in  Fig.  3  is  a  phase  contrast  image  of  a  contrast-detail  phantom 
acquired  with  120  kVp  and  a  magnification  factor  of  2.  In  Fig.  3  the  phase  contrast  is  exhibited 
as  the  highlighted  edges  of  the  test  disks.  Although  phase  contrast  generally  decreases  with 
increasing  tube  voltages,  but  this  image  demonstrates  that  the  phase  contrast  improves  image 
qualities  even  for  such  a  high  tube-voltage  of  120  kVp.  The  observer  study  on  acquired  phase 
contrast  projections  with  a  range  of  magnification  factors  from  M  =2  to  M  =3  have  validated  the 
designs  made  in  the  first  period  for  the  optimal  geometric  configurations,  which  balances  the 
conflict  requirements  for  the  x-ray  spatial  coherence  and  generating  large  diffraction  fringes  in 
images  [Wu  et  al.,  2004].  Fig.  4  shows  two  images  of  a  ACR  MAP  phantom  acquired  with  the 
prototype  system.  The  Fig.  4-(left)  is  a  phase-sensitive  projection  acquired  at  100  kVp  and  a 
magnification  factor  of  2.  For  comparison,  the  Fig.  4-(  right)  shows  the  contact  radiograph  of  the 
phantom  acquired  with  the  same  kVp  and  the  same  phantom  dose.  Both  images  were  acquired 
without  any  anti-scatter  grid.  The  loss  of  the  attenuation  contrast  associated  with  the  high  tube 
voltage  is  responsible  for  the  poor  contrast  shown  in  Fig.  4-(right),  though  the  x-ray  scatter 
contributed  to  the  image  contrast  degradation  as  well.  In  contrast,  the  phase-sensitive 
projection  image  in  Fig.  4-(left)  scores  much  better  in  the  visibilities  of  fibers,  specks  and 
masses  in  the  phantom. 

In  order  to  fully  exhibit  the  phase  contrast  of  phantom  components  and  provide  their 
projective  electron  densities,  we  conducted  a  phase  retrieval  study  to  retrieve  phantom  phase 
maps  by  using  our  novel  phase  retrieval  method.  As  an  example,  shown  in  Fig.  5  is  a  phase 
map  of  an  ACR  phantom  retrieved  by  using  our  duality-based  method  from  a  single  phase- 
sensitive  projection  acquired  with  the  prototype  system.  The  tube  voltage  used  was  120  kVp 
with  an  estimated  average  glandular  dose  of  98.9  mrad.  As  is  shown  in  Fig.  5  the  phantom 
phase  map  scores  4  fibers,  3  specks  group,  and  3.5  masses  for  the  feature-visibilities. 
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Fig.  4.  Two  images  of  a  ACR  accreditation  phantom  acquired  with  the  100  kVp  and 
the  same  dose  to  phantom,  (left)  Phase  contrast  image  of  the  phantom  with  M  =  2. 
(right)  Contact  radiograph  of  the  phantom. 


While  this  glandular 
dose  is  quite  low 
compared  to  the 
phantom  dose 
levels  of  1 50-250 
mrad  commonly 
used  in 
conventional 
mammography,  we 
expect  much  more 
dose  saving  can  be 
achieved  in  future 
experiments  with 
an  improved  control 
of  the  tube-voltage 
ramping  effects  and 
with  a  quantum- 
efficient  high 
resolution  detector 

better  matched  with  the  high  tube-voltages  employed  in  the  project.  Currently  we  are  working 
on  these  planned  improvements.  Moreover,  Fig.  5  provides  as  well  a  map  of  x-ray  phase  shifts 
for  different  phantom  features.  For  example,  from  Fig.  5  we  found  that  the  phase  shift  for  the 
acrylic  frame  of  the  phantom  is  -1 141  radians,  the  phase  shift  from  the  base  of  wax-insert 
(including  the  acrylic  substrate)  is  -1096  radians,  and  the  phase  shift  from  acrylic  body  and  the 
first  mass  is  -1106  radians.  Note  that  x-ray  phase  shifts  are  given  by  the  formula 

(p^-lr^  J  p^{s)dsBS,  is  explained  in  Section  (A),  hence  x-ray  phase  shifts  should  be  of  negative 

ray 

values.  Hence  increasing  electron  densities  or  increasing  ray  path  lengths  in  tissues  will 
generate  more  negative  phase  shift  values.  The  three  retrieved  phase  values  mentioned  above 

are  obviously  consistent  with  the  notion  of  the  phase- 
shifts.  We  have  checked  the  accuracies  of  these 
retrieved  phase  values  against  the  phase-shift 
estimates  for  60  keV  x-rays  using  the  ray  tracing 
modeling  for  a  digital  model  of  the  ACR  MAP  phantom 
that  we  built,  under  the  assumption  that  the  x-rays  in  the 
experiment  have  an  average  energy  of  60  keV.  The 
ray-tracing  based  estimates  give  the  phase  values  of  - 
923  radians,  -879  radians,  -887  radians  for  the  acrylic 
frame,  wax-insert  base  and  the  first  mass,  respectively. 
The  retrieved  phase  values  differ  by  23%  to  25%  from 
these  ray-tracing  based  estimates.  This  preliminary 
comparison  shows  that  the  accuracies  of  these 
retrieved  phase  values  are  reasonable  good,  especially 
considering  the  polychromatic  nature  of  x-rays  used  in 
the  experiments,  and  the  possible  material  differences 
between  our  digital  ACR  phantom  model  and  the  real 
ACR  MAP  phantom  used  in  experiments. 

In  summary,  in  this  period  we  conducted 
phantom  imaging  studies  with  100  kVp-140  kVp  tube 
voltages  and  retrieved  the  phantom  phase  maps  for  fully 
exhibiting  the  phase  contrast  and  providing  quantitative 


Fig.  5.  A  phase  map  of  an  ACR 
accreditation  phantom  retrieved  by  using 
the  duality-based  method  from  a  single 
phase-sensitive  projection.  The  projection 
was  acquired  with  120  kVp  and  an 
estimated  average  glandular  dose  of  98.9 
mrad. 
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information  on  the  phantom  electron  densities.  The  study  shows  that  the  high  energy  x-ray 
phase-sensitive  projections  combined  with  the  duality-based  phase  retrieval  method 
demonstrated  the  potential  of  enabling  significant  radiation  reduction  and  quantification  of  tissue 
projective  electron  densities  in  future  breast  imaging  applications. 


Key  Research  Accomplishments 


•  The  robustness  of  the  phase  retrieval  methods  is  of  critical  importance  for  limiting  and 
reducing  radiation  doses  involved  in  x-ray  phase  imaging.  In  this  period  we  have 
compared  the  AP-based  phase  retrieval  method  developed  in  the  last  period  to  the 
conventional  TIE-based  method  by  analyzing  the  phase  maps  retrieved  from  the 
experimental  images  of  a  phantom.  We  showed  that  in  these  experiments  our  AP-based 
method  is  to  be  superior  to  the  TIE-based  method  for  the  robustness  in  performing  the 
phase  retrievals.  This  finding  is  of  significance  for  developing  future  clinical  applications 
of  the  x-ray  phase  imaging. 

•  Extending  the  research  work  on  the  phase  retrieval  methods,  as  a  by-product  of  this 
project,  we  found  a  general  formula  to  quantitatively  relate  the  apparent  linear 
attenuation  coefficients  to  the  sample's  linear  attenuation  coefficients  and  refractive 
indices  in  the  phase  contrast  tomography.  This  formula  overcomes  the  gross  inaccuracy 
of  the  existing  formula  in  the  literature  in  analyzing  high-resolution  phase  contrast 
tomography,  and  it  will  be  useful  for  correctly  interpreting  and  quantifying  the  apparent 
linear  attenuation  coefficients  in  x-ray  phase  contrast  tomography  in  the  biomedical 
applications. 

•  In  this  period  we  conducted  phantom  imaging  studies  with  100  kVp  -140  kVp  tube 
voltages  as  is  planned  for  developing  low-dose  and  quantitative  x-ray  phase  imaging 
techniques.  In  the  phantom  imaging  we  validated  the  optimal  geometric  configurations 
as  determined  in  the  first  period  based  on  the  phase  contrast  visibility  analysis.  We 
performed  the  Monte  Carlo  simulations  as  well  for  accurately  estimating  the  average 
glandular  doses  involved  in  the  high-voltage  imaging.  In  order  to  fully  exhibit  the  phase 
contrast  of  phantom  components  and  provide  their  projective  electron  densities,  we 
performed  a  phase  retrieval  study  to  retrieve  phantom  phase  maps  by  using  a  novel 
phase  retrieval  method  that  we  developed.  The  results  of  the  study  show  that  the  phase 
contrast  exhibited  in  the  retrieved  phase  maps  compensates  for  the  attenuation-contrast 
loss  associated  with  using  high  tube-voltages,  and  with  these  techniques  the  radiation 
doses  are  reduced  and  quantitative  information  on  phantom  projective  electron  densities 
is  obtained.  These  studies  demonstrate  that  the  techniques  developed  here  has  good 
potential  of  enabling  significant  radiation  dose  reduction  and  quantitative  imaging  with 
tissue  electron  density  information  in  future  breast  imaging  applications. 


•  Reportable  Outcomes 


In  this  project  period  we  have  published  the  research  results  in  peer-reviewed  journals 
and  presented  the  work  in  national  conferences,  as  listed  in  the  following. 

Peer-Reviewed  Journal  Article: 
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A.  Yan,  X.  Wu,  and  H.  Liu,  Robustness  of  phase  retrieval  methods  in  x-ray  phase 
contrast  imaging:  A  comparison,  Medical  Physics  38:  5073-5080  (2011) 

A.  Yan,  X.  Wu,  Apparent  Linear  Attenuation  Coefficients  in  Phase  Contrast  X-Ray 
Tomography,  Nuclear  Instruments  and  Methods  in  Physics  Research,  B  269  1841-1843 
(2011) 


Published  Abstract  and  Conference  Presentation: 

A.  Yan,  X.  Wu,  M.  A.  Yan,  X.  Wu,  and  H.  Liu:  Development  of  a  prototype  phase- 
sensitive  x-ray  breast  imaging  system”,  DoD  Breast  Cancer  Research  Program  Era  of 
Hope  2011  Conference,  Aug.  3-Aug.  5,  2011,  Orlando,  Florida. 

H.  Liu  H,  X.  Wu:  Development  and  characterization  of  phase  sensitive  x-  ray  imaging 
systems.  The  53rd  Annual  Meeting  of  the  American  Association  of  Physicists  in 
Medicine,  July  31-Aug.  4,  201 1 ,  Vancouver,  Canada. 


Conclusion 


With  the  grant  support  by  USAMRMC  we  have  successfully  conducted  studies  on  x-ray 
phase  imaging  in  the  third  period  of  this  project.  As  the  phase  retrieval  is  a  crucially  important 
task  of  our  project,  in  this  period  we  have  examined  the  robustness  of  phase  retrieval  methods 
against  the  projection  image  noises.  In  this  period  we  compared  the  robustness  of  two  phase- 
retrieval  methods,  as  the  common  TIE-based  method  vs.  the  AP-based  method  that  we 
developed,  by  analyzing  the  phase  maps  retrieved  from  the  experimental  images  of  a  phantom. 
The  comparison  shows  that  the  conventional  TIE-based  phase  retrieval  method  is  unstable 
against  the  noise  in  the  wrap’s  projection  images,  while  the  AP-based  method  is  shown  in  the 
experiment  to  be  superior  to  the  TIE-based  method  for  the  robustness  in  performing  the  phase 
retrieval.  This  finding  is  of  importance  for  developing  the  phase  imaging  techniques. 

In  this  period  we  made  significant  progress  in  developing  low-dose  and  quantitative  x- 
ray  phase  imaging.  We  have  conducted  phantom  imaging  studies  with  100  kVp  -140  kVp  tube 
voltages  and  retrieved  the  phantom  phase  maps  for  fully  exhibiting  the  phase  contrast  and 
providing  quantitative  information  on  the  phantom  electron  densities.  The  results  of  the  study 
show  that  the  phase  contrast  exhibited  in  the  retrieved  phase  maps  compensates  for  the 
attenuation-contrast  loss  associated  with  using  high  tube-voltages,  and  with  these  techniques 
the  radiation  doses  are  reduced  and  quantitative  information  on  phantom  projective  electron 
densities  is  obtained.  These  techniques  developed  in  this  period  has  good  potential  of  enabling 
significant  radiation  dose  reduction  and  quantitative  imaging  with  tissue  electron  density 
information  in  future  breast  imaging  applications. 

In  the  next  no-cost  extension  period,  we  will  continue  our  research  in  developing  low- 
dose  and  quantitative  x-ray  phase  imaging  with  high  energy  x-rays.  We  plan  to  further  optimize 
the  design  trade-off  in  the  detector  spatial  resolutions  and  quantum  efficiencies  for  achieving 
more  large  radiation  dose  reductions.  We  will  refine  our  phase  retrieval  method  for  higher 
retrieval  accuracies  by  better  incorporating  the  x-ray  spectral  effects.  These  studies  in  next 
period  will  facilitate  the  development  of  the  low  dose  quantitative  phase  imaging  techniques  for 
future  breast  imaging  applications. 
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Purpose:  The  robustness  of  the  phase  retrieval  methods  is  of  critical  importance  for  limiting  and 
reducing  radiation  doses  involved  in  x-ray  phase  contrast  imaging.  This  work  is  to  compare  the 
robustness  of  two  phase  retrieval  methods  by  analyzing  the  phase  maps  retrieved  from  the  experi¬ 
mental  images  of  a  phantom. 

Methods:  Two  phase  retrieval  methods  were  compared.  One  method  is  based  on  the  transport  of  in¬ 
tensity  equation  (TIE)  for  phase  contrast  projections,  and  the  TIE-based  method  is  the  most  com¬ 
monly  used  method  for  phase  retrieval  in  the  literature.  The  other  is  the  recently  developed 
attenuation-partition  based  (AP-based)  phase  retrieval  method.  The  authors  applied  these  two  meth¬ 
ods  to  experimental  projection  images  of  an  air-bubble  wrap  phantom  for  retrieving  the  phase  map 
of  the  bubble  wrap.  The  retrieved  phase  maps  obtained  by  using  the  two  methods  are  compared. 

Results:  In  the  wrap’s  phase  map  retrieved  by  using  the  TIE-based  method,  no  bubble  is  recogniz¬ 
able,  hence,  this  method  failed  completely  for  phase  retrieval  from  these  bubble  wrap  images.  Even 
with  the  help  of  the  Tikhonov  regularization,  the  bubbles  are  still  hardly  visible  and  buried  in  the 
cluttered  background  in  the  retrieved  phase  map.  The  retrieved  phase  values  with  this  method  are 
grossly  erroneous.  In  contrast,  in  the  wrap’s  phase  map  retrieved  by  using  the  AP-based  method, 
the  bubbles  are  clearly  recovered.  The  retrieved  phase  values  with  the  AP-based  method  are  reason¬ 
ably  close  to  the  estimate  based  on  the  thickness-based  measurement.  The  authors  traced  these  stark 
performance  differences  of  the  two  methods  to  their  different  techniques  employed  to  deal  with  the 
singularity  problem  involved  in  the  phase  retrievals. 

Conclusions:  This  comparison  shows  that  the  conventional  TIE-based  phase  retrieval  method, 
regardless  if  Tikhonov  regularization  is  used  or  not,  is  unstable  against  the  noise  in  the  wrap’s 
projection  images,  while  the  AP-based  phase  retrieval  method  is  shown  in  these  experiments 
to  be  superior  to  the  TIE-based  method  for  the  robustness  in  performing  the  phase  retrieval. 

©2011  American  Association  of  Physicists  in  Medicine.  [DOI:  10.1118/1.3618731] 

Key  words:  x-ray  imaging,  phase  contrast,  phase  retrieval 


I.  INTRODUCTION 

In  recent  years,  the  in-line  phase-contrast  x-ray  imaging  has 
attracted  intensive  research  efforts.  The  in-line  phase-con¬ 
trast  x-ray  imaging  is  a  technique  using  the  free  space  dif¬ 
fraction  of  phase-shifted  x-rays  to  form  the  interference 
fringes  at  tissues’  boundaries  and  interfaces  in  images.'”*^  In 
fact,  for  over  100  years,  the  tissue  attenuation  differences 
have  been  the  sole  contrast  mechanism  for  medical  x-ray 
imaging.  However,  when  x-rays  traverse  the  body  parts,  as  a 
wave  x-rays  undergo  phase  shifts  as  well.  The  amount  of  the 
phase  shift  along  an  exit  ray  is  determined  by  the  tissue  re¬ 
fractive  indices  along  the  ray.  Note  that  x-ray  refractive 
index  n  is  a  complex  number  and  equal  to  «  =  1  —  (5  -f  //?, 
where  5  is  the  refractive  index  decrement  and  responsible 
for  x-ray  phase  shift,  while  fl  is  the  imaginary  part  of  the  re¬ 
fractive  index  and  responsible  for  x-ray  attenuation.  The 
amount  of  x-ray  phase  shift  along  an  exit  ray  is  given  by 
(^  =  —{2n/ a)  ^  6ds,  where  X  is  the  x-ray  wavelength,  and 


the  integral  is  over  the  ray  path.^’^  In  other  words,  the  phase 
shift  is  equal  to  the  projection  of  refractive  index  decrements 
scaled  by  a  factor  {In/ X).  On  the  other  hand,  x-ray  attenua¬ 
tion  along  an  exit  ray  is  determined  by  the  projection  of  the 
tissue  linear  attenuation  coefficients  along  the  ray,  which  is 
equal  to  {An/ X)  J  flds.^’^  We  have  estimated  S  and  /(  values 
for  the  biological  tissues  and  found  that  the  tissue  5  values 
(10~®-10~*^)  are  about  1000  times  greater  than  their  /?  values 
(10“^-10“^')  for  10-150  keV  x-rays.*  Hence,  the  differen¬ 
ces  in  x-ray  phase  shifts  between  different  tissues  are  about 
1000  times  greater  than  their  differences  in  the  projected  lin¬ 
ear  attenuation  coefficients.  Therefore,  the  phase-contrast 
imaging  techniques  may  greatly  increase  the  lesion-detection 
sensitivity  for  x-ray  imaging.  The  settings  for  the  inline 
phase-sensitive  imaging  are  similar  to  that  of  conventional 
x-ray  imaging,  provided  a  source  with  very  small  focal  spot 
and  a  sufficiently  large  object-detector  distance  are 
required.^’^’*  In  the  inline  imaging  setting,  x-rays  undergo 
phase  shifts  as  traversing  the  imaged  object,  and  then  diffract 
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freely  over  a  sufficiently  large  distance  before  reaching  the 
detector.  In  this  way,  the  tissues’  phase  contrast  manifests  as 
the  dark-bright  diffraction  fringes  at  tissues’  boundaries  and 
interfaces  in  the  measured  images.^’^’*  Hence,  the  inline 
phase  contrast  imaging  has  good  potential  of  greatly  enhanc¬ 
ing  the  detection  sensitivity  and  reducing  radiation  doses 
involved  in  the  imaging. 

However,  as  the  interfaces  and  boundaries  of  the  different 
tissue  compartments  are  greatly  accentuated  in  a  phase-con¬ 
trast  image,  the  bulk  phase  contrast  in  a  given  tissue  com¬ 
partment,  where  the  phase  shifts  may  vary  slowly,  could  get 
lost.  This  is  because  the  phase  contrast  is  proportional  to  the 
Laplacian  and  gradient  differentials  of  the  phase  shifts,  as  is 
shown  in  Sec.  II  below.  Moreover,  the  information  about  x- 
ray  phase  shifts  are  valuable  for  tissue  characterizations, 
since  x-ray  phase  shift  along  a  ray  is  proportional  to  the  pro¬ 
jected  tissue  electron  density,  that  is,  (p  =  —h-f,  J  p^ds, 
where  ;•(,  is  the  classic  electron  radius  and  denotes  the  tis¬ 
sue  electron  density  and  the  integral  is  over  the  ray  path.*’^ 
In  order  to  fully  exhibit  tissue  phase  contrast  and  reconstruct 
tissue  projected  electron  densities  for  quantitative  tissue 
characterizations,  one  needs  to  extract  the  tissue  phase  shifts 
from  the  mixed  contrast  exhibited  in  a  phase-sensitive  pro¬ 
jection.  The  procedure  of  retrieving  the  phase-shift  map  of 
an  object  from  its  phase-sensitive  projections  is  called  the 
phase  retrieval.  A  retrieved  phase  map  is  able  to  provide 
a  quantitative  map  of  the  object’s  projected  electron  den¬ 
sities,  which  could  be  used  for  quantitative  tissue  character¬ 
izations.^^”*^  Moreover,  performing  phase  retrieval  is 
necessary  for  reconstructing  volumetric  3-D  maps  of  tissue 
attenuation  coefficients  and  refractive  indices,  respec¬ 
tively,*"*’*^  and  for  eliminating  the  phase-contrast  caused 
artifacts  in  the  volumetric  3-D  images.*® 

Examining  the  robustness  of  phase  retrieval  methods 
against  the  projection  noise,  such  as  the  x-ray  photon  quan¬ 
tum  noise  and  others  involved  in  the  projection  acquisitions, 
is  an  important  task  in  phase  contrast  imaging.  If  a  phase  re¬ 
trieval  method  is  not  robust  against  the  noise,  the  method 
could  become  unstable  and  fail  in  phase  retrievals  if  the 
acquired  projection  images  had  be  contaminated  with  sub¬ 
stantial  noise.  On  the  other  hand,  suppressing  x-ray  quantum 
noise  may  require  using  high  radiation  doses  in  the  image 
acquisitions.  For  clinical  applications,  it  is  imperative  to 
limit  and  reduce  radiation  doses  involved,  hence,  it  is  critical 
to  develop  robust  phase  retrieval  methods  for  future  clinical 
applications  of  phase  contrast  imaging.  The  purpose  of  this 
work  is  to  compare  the  robustness  of  two  phase  retrieval 
methods  by  means  of  analyzing  the  phase  maps  retrieved 
from  the  experimental  projection  images  of  an  air-bubble 
wrap  phantom.  We  will  first  introduce  the  phase  retrieval 
method  that  is  based  on  the  transport  of  intensity  equation 
(TIE),  which  describes  how  the  phase  contrast  is  encoded  in 
the  projection  images.***’**  This  TIE-based  method  is  the 
most  commonly  used  method  for  phase  retrieval  in  the  litera¬ 
ture.  We  will  then  introduce  a  recently  developed  phase  re¬ 
trieval  method,  namely  the  attenuation-partition  based  (AP- 
based)  method.*^  In  order  to  examine  the  robustness  of  these 
two  methods  against  noises,  we  applied  these  two  phase  re¬ 


trieval  methods  to  experimental  phantom  images,  namely  a 
radiograph  and  a  phase  contrast  image  of  the  air-bubble 
wrap  phantom,  for  the  phase  retrievals.  We  will  analyze  the 
retrieved  phase  maps  obtained  by  using  the  two  methods  and 
compare  their  performances  in  the  robustness  against  noises 
in  the  phantom  images.  In  Sec.  IV,  we  explain  that  these  per¬ 
formance  differences  of  the  two  phase  retrieval  methods  are 
rooted  in  their  different  techniques  employed  to  deal  with 
the  singularity  problem  involved  in  the  phase  retrievals.  In 
addition,  other  factors  that  may  affect  the  phase  retrieval  per¬ 
formance  will  be  briefiy  discussed  as  well. 

II.  PHASE  RETRIEVAL  METHODS 

In  phase  contrast  images,  the  attenuation  contrast  and 
phase  contrast  are  mixed  together.  In  order  to  retrieve  phase 
maps  from  phase  contrast  projection  images,  one  should 
understand  how  the  phase  contrast  is  encoded  in  the  projec¬ 
tion  images.  This  understanding  can  be  gained  from  the  x- 
ray  propagation  equation  such  as  the  TIE.***’**  This  equation 
can  be  derived  either  from  the  x-ray  Fresnel-diffraction 
equation  or  the  Wigner  distribution  based  phase-space  for¬ 
malism.^”**  If  we  denote  the  x-ray  transmission  image  of  an 
object,  or  equivalently  its  x-ray  attenuation  map  by  A^(r) 
and  the  x-ray  phase-shift  map  of  the  object  by  then  the 
detected  x-ray  intensities  are  given  by**’**’**' 

=  w  V  •  (A2(,-^/M)V^(f),/M))| 

(1) 

where  X  is  the  average  x-ray  wavelength  of  the  polychro¬ 
matic  x-rays.  In  Eq.  (1),  is  the  entrance  x-ray  intensity, 
Ri  is  the  source-object  distance  set  in  the  projection,  R2  is 
the  object-detector  distance,  M  =  {Ri  +  R2)/Ri,  the  magni¬ 
fication  factor  in  the  projection,  and  Fq  denotes  the  position 
vector  in  the  detector  plane.  In  Eq.  (1),  the  operator  V 
denotes  the  two-dimensional  transverse  gradient  differential 
operator.  For  the  purpose  of  this  work,  we  assume  that  the 
x-ray  source’s  focus  spot  is  ideally  point-like  and  the  detec¬ 
tor  employed  is  an  ideal  detector.  Obviously,  the  detected 
x-ray  intensity  /(fp)  is  determined  not  only  by  attenuation 
map  Ag{r)  in  the  projection,  but  also  by  the  encoded  phase 
contrast,  that  is,  by  the  transverse  Laplacian  and  gradient 
differentials  of  the  phase-shift  map  (p{r)  in  the  projection. 
Since  the  phase  contrast  and  attenuation  contrast  are  mixed 
together  in  a  phase  contrast  image,  as  is  shown  by  Eq.  (1), 
hence  measuring  a  single  phase-contrast  image  is  not 
enough  for  retrieving  the  phase-shift  map  (p{r)  in  the  pro¬ 
jection.  Therefore,  in  general,  at  least  two  projection  images 
are  needed  for  retrieving  the  phase-shift  map  (p{r)  of  the 
object. 

The  most  commonly  used  phase  retrieval  method  in  the  lit¬ 
erature  is  the  TIE-based  method.  In  this  method,  two  acquired 
projection  images  of  an  object  are  used  for  retrieving  the 
phase-shift  map  of  the  object:  one  is  a  contact  radiograph 
Ag{r)  of  the  object,  the  other  is  a  phase  contrast  projection 
image  /(fp)  of  the  object.  One  then  is  able  to  retrieve  the 
phase-shift  map  (p(f)  by  solving  Eq.  (1)  for  (p{r)  as  follows*^: 
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(p{F)  =  -  (27rM/;j?2)V^2 


X 


(2) 


In  above  equation,  the  operator  denotes  the  inverse  of 
the  Laplacian  differential  operator  =  {d^/dx^  +  d^/dy^). 
The  inverse  Laplacian  operator  is  a  pseudodifferential 
operator.  The  action  on  a  function  g{r)  of  a  pseudodifferential 
operator  such  as  is  dehned  as 


exp(/27t(T-  f)  -f)  ■  (-1/(471^/^)) 
X  cf-f. 


(3) 


In  Eq.  (3),  s  and/  denote  the  integral  variables  in  the  coordi¬ 
nate  space  and  frequency  space,  respectively.  Hence,  one 
can  compute  the  phase  map  (p{f)  by  using  Eq.  (2)  with  the 
help  of  Eq.  (3).  The  TIE-hased  method  is  effective  and  fast 
for  phase  retrievals  in  the  cases  with  strong  phase  contrast 
effects  and  low  noise  levels  as  is  shown  in  many  cases  dis¬ 
cussed  in  the  literature. However,  the  Achilles  heel  of 
the  TIE-based  method  lies  at  the  inverse  Laplacian  operator 
involved  in  Eq.  (2).  The  operator  adds  a  zero-fre¬ 
quency  singularity  to  Eq.  (2),  as  is  suggested  by  Eq.  (3).  As 
we  will  see  below,  the  singularity  may  amplify  the  noise  ran¬ 
domly  embedded  in  the  acquired  images  and  result  in  insta¬ 
bility  in  phase  retrievals. 

In  order  to  get  rid  of  the  singular  pseudodifferential  oper¬ 
ator  involved  in  the  TIE-based  phase  retrievals,  we 
recently  developed  a  novel  phase  retrieval  method:  the  AP 
based  method  or  the  AP-based  method  for  short.*^’'*^  Our 
idea  in  this  method  is  to  utilize  the  correlation  between  the 
x-ray  phase  shift  and  its  attenuation  to  eliminate  any  singu¬ 
larity  involved  in  the  phase  retrievals.  As  is  well  known,  in 
the  diagnostic  x-ray  imaging,  x-ray  attenuation  arises 
from  three  x-ray-matter  interactions:  the  photoelectric 
absorption,  the  coherent  scattering,  and  the  incoherent  scat¬ 
tering.  Correspondingly,  we  can  partition  the  x-ray  attenua¬ 
tion  into  two  parts:  the  Compton  scattering-caused 
attenuation  and  the  attenuation  caused  hy  photoelectric 
absorption  and  coherent  scattering,  which  is  denoted  by 
A^ecohi^-  That  is,  we  can  partition  the  overall  x-ray  attenua¬ 
tion  A^  into  a  product  of  two  parts 


(4) 


The  reason  of  factoring  out  A^^  is  that  both  of  A^^  and  the 
phase-shift  cp  are  determined  by  the  tissue’s  electron  density 


Akn{i^  =  exp 


(p(r)  =  -AfeP.^pir)  (5) 


where  =  2.818  x  10^^^  m  is  the  classical  electron  radius, 
Pgp,  the  projected  electron  density  along  the  ray  path,  that  is, 
Pepif)  —  Okn  denotes  the  average  Comp¬ 

ton  scattering  cross-section  for  polychromatic  x-rays.  Note 
that  Compton  scattering  cross-section  is  given  by  the  Klein- 
Nishina  total  cross-section*^ 


I  T  L 


2(1  +  ;?) 


1  +  2r\ 


-log(l  +2t]) 
'7 


-f  — log(l-|-27;) 
2ri 


1  +  3?7  ] 

(1  +  2^)"/’ 


(6) 


where  t]  =  E  /511  keV  and  it  slowly  varies  with  x-ray 
energy  E  for  E^  511  keV.  We  observed  that  the  extent  of 
the  correlation  between  phase  and  attenuation  depends  on 
the  x-ray  photon  energies  and  the  tissues  elemental  composi¬ 
tions.  For  example,  for  light  elements  with  atomic  numbers 
Z  <  10  and  x-rays  of  60  keV  or  higher,  the  x-ray-matter 
interactions  are  dominated  by  the  Compton  scattering, 
hence,  both  the  tissue  attenuation  and  phase  shift  are  all 
determined  by  tissues’  electron  density  distributions.  We  call 
this  relationship  between  phase  shift  and  attenuation  the 
phase-attenuation  duality,  and  we  dehne  the  so-called  duality 
transform  as^** 

=  (1  -  {I^R2re/2TiMaKN)y^y'-(^^^^y 

(7) 

For  the  cases  where  the  x-ray-matter  interactions  are  domi¬ 
nated  by  the  Compton  scattering,  we  proved  that  the  attenua¬ 
tion  A|^  can  be  found  from  a  single  phase  contrast  image 
7(Td)  by  means  of  the  duality  transform  such  that 
=  'D(/(fo)).^**  With  the  help  of  the  phase-attenuation  duality, 
the  phase-shift  map  of  the  object  can  be  simply  retrieved  as 

(p{r)  =  ^lnAi^(F)  =  ^In  35(/(Fo)).  (8) 


We  call  this  formula  as  the  phase-attenuation  duality-based 
phase  retrieval  formula.^**  One  significant  advantage  of  this 
duality-based  phase  retrieval  formula  is  that  the  pseudodif¬ 
ferential  operator  (l  —  (^)?R2i'e/2nMd'KN)'’^^)  involved  in 
Eqs.  (7)  and  (8)  is  free  of  any  singularity.  Therefore,  the 
phase  retrieval  method  based  on  the  phase-attenuation  dual¬ 
ity  is  stable  and  robust  against  the  noise  in  images.^** 

Unfortunately,  the  phase-attenuation  duality  does  not  hold 
in  some  other  cases  such  as  imaging  with  low  energy  x-rays  or 
imaging  bones  and  calcifications  that  contain  substantial 
amount  of  heavy  elements.  In  those  cases,  the  retrieved  phase 
map  using  Eq.  (8)  may  provide  merely  coarse  estimates  of  the 
real  phase  maps  of  the  imaged  objects.  In  order  to  overcome 
this  limitation  of  the  duality-based  phase  retrieval  formula  Eq. 
(8),  recently  we  proposed  to  correct  the  errors  iteratively 
through  repeated  comparisons  of  the  computed  estimates  of 
the  x-ray  Fresnel-diffraction  intensities  against  the  acquired 
projection  intensities.  To  implement  this  strategy,  we  devel¬ 
oped  the  attenuation-partition  based  iterative  algorithm,  whose 
flow  chart  is  shown  in  Fig.  1.*^  While  interested  readers  can 
find  the  mathematical  proofs  of  the  algorithm  in  Ref.  13,  here, 
we  give  a  brief  explanation  of  this  algorithm  flow  chart  in 
Fig.  1.  In  this  flow  chart,  A^  denotes  the  attenuation  map  of 
the  imaged  object  measured  from  its  contact  radiograph  and  I 
is  the  acquired  phase  contrast  image  of  the  object.  In  the  flow 
chart,  Ikn  denotes  the  simulated  phase  contrast  image  formed 
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Fig.  1 .  Flow  chart  of  the  AP  based  iterative  algorithm. 

by  using  the  estimated  object’s  electron  densities  and  by 
turning  off  the  attenuation  processes  associated  with  x-ray 
photoelectric  absorption  and  coherent  scattering.  In  Fig.  1,  two 
transforms  are  implemented  in  the  iterations,  one  is  the  Fres- 
nel-diffraction  transform  which  implements  x-ray 

Fresnel-diffraction,^'  and  the  other  is  the  duality  transform 
as  is  defined  in  Eq.  (7).  According  to  this  flow  chart  (Fig.  1)  of 
the  AP-based  phase  retrieval  method,  one  first  applies  the 
duality  transform  T)  to  the  acquired  phase  contrast  image  I 
to  obtain  an  estimate  of  and  the  phase  map  cp  of  the 
object.  One  then  compares  A|^  against  the  measured  A^ 
and  computes  the  error  5 A  =  Ay^jv  —  Aq.  Applying  the  Fres- 
nel-diffraction  transform  gre  to  the  fictitious  transmission 
function  <5A  •  exp{i(p),  one  computes  the  Ikn  -correction 
51  =  1 •  e\p{i(p))  I  ,  and  a  new  estimate  of  I^n  is  deter¬ 
mined  as  /jHv  =  (V^  +  With  this  new  estimate  of 

one  can  start  a  new  round  of  the  iteration  as  is  indicated  in  the 
flow  chart  Fig.  1.  The  iteration  converges  when  the  does 
not  change  substantially  with  further  iteration  steps. 


III.  COMPARISON  OF  THE  TWO  PHASE  RETRIEVAL 
METHODS 

In  order  to  compare  the  robustness  of  above  two  phase  re¬ 
trieval  methods,  namely,  the  TIE-based  method  and  the  AP- 
based  method,  we  applied  the  two  methods  for  retrieving  the 
phase  map  of  a  piece  of  air-bubble  wrap.  In  the  bubble  wrap, 
the  air  packets  are  locked  between  two  thin  layers  of  low-den¬ 
sity  polyethylene  hlms  forming  the  air  bubbles.  Shown  in  Fig. 
2(a)  is  a  contact  radiograph  of  the  wrap.  The  x-ray  source 
used  was  a  micro-focus  x-ray  tube  (L8 121-02,  Hamamatsu) 
with  a  focal  spot  size  of  7  fim.  The  source  has  a  tungsten  tar¬ 
get  without  added  filtration  and  was  operating  at  40  kVp  and 
0.2  mA  for  a  30  s  exposure.  The  average  x-ray  photon  energy 
was  estimated  at  14.6  keV  for  the  40  kVp  x-rays.  In  this  ac¬ 
quisition,  the  source-detector  distance  (SID)  was  set  to  1.78 
m.  The  imaging  detector  employed  was  an  aSe-based  flat- 
panel  detector  (DirectRay,  Hologic)  with  a  pixel  pitch  of  140 
^m.  Since  the  thin  polyethylene  films  and  locked  air  bubbles 
in  fhe  wrap  generate  little  differences  in  x-ray  attenuation, 
hence,  the  wrap’s  radiograph  in  Fig.  2(a)  exhibits  only  little 
image  contrast,  and  the  x-ray  quantum  noise  in  the  acquisition 
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Fig.  2.  (a)  Contact  radiograph  of  the  bubble  wrap  acquired  at  40  kVp  and 
with  a  SID  =  1 .75  m.  (b)  Intensity  profile  along  the  marked  line  on  (a). 


masked  the  wrap’s  details  such  that  the  rims  of  many  bubbles 
are  not  visible  on  the  image.  Figure  2(b)  shows  a  profile  of 
image  intensifies  along  fhe  marked  dash  line  on  Fig.  2(a)  and 
fhe  prohle  reveals  how  the  noisy  intensity  variations  mask  the 
low  contrast  rims  of  the  bubbles.  In  contrast.  Fig.  3(a)  is  a 
phase  contrast  image  of  the  bubble  wrap.  This  phase  contrast 
image  was  acquired  with  a  sample-detector  distance 
7?2  =  1.15  m  and  otherwise  the  same  technique  settings  as 
that  used  for  the  radiograph,  that  is,  with  the  same  focal  spot 
size,  same  tube  voltage,  same  tube  current  and  exposure  time, 
same  SID,  and  the  same  detector.  However,  in  this  projection, 
the  phase-shifted  x-rays  were  allowed  to  freely  diffract  over 
1.1 5-m  long  distance  on  their  way  to  the  detector  and  to  form 
the  interference  fringes  at  the  bubble  boundaries  on  the  image. 
As  is  shown  in  Fig.  3(a),  the  phase  contrast  depicts  not  only 
the  bubble  rims,  but  also  the  dents  inside  the  bubble  domes. 
The  enhancement  is  resulted  from  the  rapid  changes  of  the 
projected  thicknesses  of  the  polyethylene  films  at  the  bubble 
rims  and  the  dents  inside  bubble  domes.  These  rapid  changes 
in  the  projected  thicknesses  of  polyethylene  generate  as  well 
rapid  changes  in  x-ray  phase  shifts,  since  the  x-ray  phase 
shifts  are  given  by  (p{r)  =  —XrePgTp{f),  where  denotes 
polyethylene  film’s  elecfron  density  and  Tp{f)  denotes  pro¬ 
jected  thickness  of  polyethylene  along  the  ray.  According  to 
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Fig.  3.  (a)  Phase  contrast  image  of  the  bubble  wrap,  which  was  acquired  at 
40  kVp  and  with  a  SID  =1.75  m  and  a  magnification  factor  M  =  2.8. 
(b)  Intensity  profile  along  the  marked  line  on  Fig.  3(a). 


Eq.  (1),  the  Laplacian  and  gradient  differentials  of  the  rapid 
phase  shifts  generated  large  x-ray  intensity  variations,  which 
manifest  as  the  bright-dark  fringes  flank  along  the  bubble  rims 
in  Fig.  3(a).  Figure  3(b)  shows  a  profile  of  image  intensities 
along  the  marked  dash  line  on  Fig.  3(a),  as  the  line  traces  the 
same  positions  as  indicated  by  the  dash  line  on  Fig.  2(a).  The 
downward  and  upward  overshooting  of  the  intensity  values  in 
the  profile  represent  the  bright  and  dark  fringes  at  the  bubble 
boundaries  in  phase  contrast  image  Fig.  3(a).  These  up-down 
bipolar  overshootings  of  intensities  are  much  larger  than  the 
background  noise,  hence,  the  noise  gets  masked  in  Fig.  3(a). 
While  the  edge-enhancement  generated  by  phase-contrast  is 
generally  useful  for  imaging  the  wrap,  however,  such  edge- 
enhancements  may  lead  interpretation  errors  in  the  characteri¬ 
zation  of  the  wrap’s  structure  and  composition.  For  example, 
in  Fig.  3(a),  the  dark  fringes  flanked  the  enhanced  bubble  rims 
may  be  confused  with  possible  structural  breaks  in  the  wrap 
sample,  although  we  knew  in  advance  that  this  wrap  sample  is 
free  of  any  damage.  In  order  to  fully  exhibit  the  wrap’s  phase 
contrast  and  correctly  characterize  the  bubble  wrap,  it  is  nec¬ 
essary  to  perform  the  phase  retrieval.  In  order  to  retrieve  the 
wrap’s  phase-shift  map  from  the  wrap’s  radiograph  [Fig.  2(a)] 
and  its  phase  contrast  image  [Fig.  3(a)],  we  first  employed 
the  TIE -based  phase  retrieval  formula  Eq.  (2).  The  retrieved 


phase  map  with  the  TIE-method  is  shown  Fig.  4(a).  Appa¬ 
rently,  the  phase  retrieval  failed  completely,  since  no  bubble 
is  recognizable  in  Fig.  4(a).  Although  the  failure  can  be 


Fig.  4.  Retrieved  phase  maps  of  the  bubble  wrap,  (a)  With  the  TIE-based 
method,  (b)  With  the  TIE-based  method  and  Tikhonov  regularization,  (c) 
Profile  of  the  retrieved  phase  values  along  the  marked  line  on  (b). 


Medical  Physics,  Vol.  38,  No.  9,  September  2011 


5078 


Yan,  Wu,  and  Liu:  Robustness  of  phase  retrieval  methods 

ascribed  to  the  high-levels  of  the  quantum  noise  presented  in 
the  radiograph  of  Fig.  2(a),  but  it  is  rooted  in  the  intrinsic 
instability  of  the  TIE-method.  As  we  pointed  out  in  Sec.  II, 
the  inverse  Laplacian  operator  plugs  in  the  zero-fre¬ 
quency  singularity  into  the  TIE-base  phase  retrieval  formula 
Eq.  (2).  The  singularity  amplified  the  noise  randomly  embed¬ 
ded  in  the  projections  and  results  in  the  failure  in  the  phase 
retrievals.  To  mend  the  instability  of  the  TIE-method,  one  can 
try  to  apply  the  regularization  techniques  to  tame  the  singular¬ 
ity-caused  problems.  A  common  regularization  technique  in 
the  literature  is  Tikhonov  regularization.^^  With  this  regulari¬ 
zation  scheme,  one  replaces  the  inverse  Laplacian  operator 
by  a  pseudodifferential  operator  V^/[(V^)^  -t-  a^],  where 
a  is  the  regularization  parameter,  which  is  roughly  propor¬ 
tional  to  the  images  noise-signal  ratios  of  the  images.^^  In 
essence,  Tikhonov  regularization  seeks  the  minimum-norm, 
least  squares  solution  of  Eq.  (2).  While  the  Tikhonov  regulari¬ 
zation  of  the  inverse  Laplacian  operator  may  tame  the  insta¬ 
bility,  it  sacrifices  the  phase  retrieval  accuracies.  The 
performance  in  phase  retrieval  of  Tikhonov  regularization  is 
highly  dependent  on  the  amounts  of  noise  presented  in  the 
images.  The  regularization  parameter  a  was  selected  through 
a  trial  by  comparing  the  phase  maps  retrieved  with  a  wide 
range  of  a-values  such  as  a  =  Am^,  5Am^,  10  Am^,  50  Am^, 
100  Am^,  200  Am^,  300  Am^,  and  500  Am^,  where  Am  is  the  fre¬ 
quency  sampling-step  used  in  the  phase  retrieval.  Figure  4(b) 
is  the  wrap  phase  map  retrieved  using  the  TIE-method  and 
Tikhonov  regularization  with  a  =  200  and  Au^  —  3.782 
X  10^^  which  was  the  a  -value  for  the  best  results  as 

determined  from  the  trial.  Apparently,  these  retrieval  results 
are  very  unsatisfactory.  In  Fig.  4(b),  the  bubble  rims  are 
hardly  visible  in  the  extremely  cluttered  background.  Figure 
4(c)  shows  the  profile  of  the  retrieved  phase  values  along  the 
marked  dash  line  on  Fig.  4(b).  As  we  will  explain  below  that 
these  phase  values  are  grossly  erroneous. 

As  a  comparison,  we  have  applied  the  AP-based  method 
to  the  wrap’s  images  in  Figs.  2(a)  and  3(a)  for  the  phase  re¬ 
trieval  as  well.  Following  the  method’s  iteration  flow  chart 
in  Fig.  1,  we  succeeded  in  retrieving  the  wrap’s  phase  map 
of  the  bubble  wrap  with  ten  iteration  steps.  Figure  5(a)  is  the 
retrieved  phase  map  of  the  wrap  with  the  AP-based  method. 

In  a  stark  contrast  to  Figs.  4(a)  and  4(b)  where  the  bubble 
rims  are  hardly  visible,  the  bubble  rims  are  prominently 
depicted  in  Fig.  5(a).  A  profile  of  retrieved  phase  values 
along  the  marked  line  on  Fig.  5(a)  is  shown  in  Fig.  5(b). 
Note  that  all  the  marked  lines  in  Figs.  3(a),  4(b),  and  5(a), 
respectively  correspond  to  the  same  line  defined  on  the  bub¬ 
ble  wrap.  In  this  way,  one  can  easily  compare  the  three  pro¬ 
files  shown  in  Figs.  3(b),  4(c),  and  5(b).  In  order  to  study  the 
quantitative  aspects  of  the  retrieved  phase  maps  of  the  wrap, 
remember  that  the  amount  of  phase  shift  along  a  ray  is  given 
by  (p{f)  =  —XrepJTeif),  here,  denotes  polyethylene  film’s 
electron  density  and  Tp{f)  denotes  projected  thickness  of 
polyethylene  along  the  ray.  Note  that  x-ray  phase  shifts 
should  be  of  negative  values,  as  x-ray  refractive  indices  of 
tissues  and  materials  are  complex  and  their  real  part  are  less 
than  one.  The  x-rays  traversed  longer  paths  in  polyethylene 
at  the  bubble  rims  and  incurred  larger  phase  shifts,  as 
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Fig.  5.  (a)  Retrieved  phase  maps  of  the  bubble  wrap  with  the  AP-based 
method,  (b)  Profile  of  the  retrieved  phase  values  along  the  marked  line  on  (a). 


compared  to  other  parts  of  the  wrap.  For  example,  the  three 
sharp  negative  peaks  in  the  profile  Fig.  5(b)  reflect  the  large 
x-ray  phase  shifts  incurred  at  the  three  rim  locations  along 
the  marked  line  in  Fig.  5(a).  In  order  to  gauge  the  accuracies 
of  the  retrieved  phase  values  profiled  in  Fig.  5(b),  one  has  to 
know  the  values  of  Tp{f),  the  projected  thicknesses  in  poly¬ 
ethylene  films.  Apparently,  Tp{f)  depends  on  the  local  curva¬ 
tures  of  bubble  domes  and  the  incident  angles  of  the  rays. 
While  it  is  hard  to  measure  Tp{r)  values  in  the  bubbles,  we 
measured  the  thicknesses  of  the  flat  bases  between  bubbles 
using  a  caliber  ruler,  and  we  found  that  the  base  thickness 
was  about  0.055  mm.  Knowing  that  the  molecular  formula 
of  the  low-density  polyethylene  film  is  (€2114)^^  and  its  mass 
density  is  0.925  g/cm  ,  we  calculated  its  electron  density  as 
3.18  X  10^^ /cm^.  Assuming  approximately  all  rays  have 
normal  incidence,  we  hnd  that  the  approximate  phase  shifts 
at  the  flat  bases  between  bubbles  is  about  —4.2  radian.  On 
the  other  hand,  according  to  the  phase  profile  in  Fig.  5(b), 
the  average  phase  shift  at  the  flat  bases  between  bubbles  is 
—  1.05  radian.  Obviously,  these  two  estimates  of  the  phase 
values  though  different,  but  are  reasonably  close.  In  compar¬ 
ison,  the  phase  profile  in  Fig.  4(c),  which  was  obtained  by 
using  the  TIE-based  method  with  Tikhonov  regularization, 
depicts  just  messy  up-down  peaks  buried  in  cluttered  and 
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noisy  background.  In  the  phase  profile,  Fig.  4(c)  only  two 
negative  peaks  at  the  rims  are  barely  recognizable,  and  the 
third  expected  negative  peak  completely  disappears.  Accord¬ 
ing  to  this  profile,  the  phase-shift  values  at  the  bases  fluctu¬ 
ate  over  a  wide  range  from  0  to  +87.5  radian  with  an 
average  of  +27.3  radian,  while  x-ray  phase  shifts  should  be 
negative.  Hence,  the  average  phase  value  of  flat  bases  esti¬ 
mated  based  on  the  profile.  Fig.  4(c)  is  in  a  gross  discrepancy 
with  the  thickness  measurement-based  estimate.  Therefore, 
the  above  comparisons  show  clearly  that  the  AP-based  phase 
retrieval  method  is  shown  to  be  superior  to  the  TIE-based 
method,  regardless  if  the  Tikhonov  regularization  is  used,  in 
performing  the  phase  retrieval  for  the  bubble  wrap. 

Before  ending  this  section,  we  want  to  demonstrate  the  use¬ 
fulness  of  performing  the  phase  retrieval  for  tissue  and  mate¬ 
rial  characterizations  by  comparing  the  bubble  wrap’s  phase 
map  in  Fig.  5(a)  to  its  phase  contrast  image  in  Fig.  3(a).  As  we 
pointed  out  earlier,  the  dark  fringes  which  are  flanking  the 
bubble  rims  in  the  phase  contrast  image  Fig.  3(a),  may  be  con¬ 
fused  for  indicating  possible  structural  breaks  in  the  bubble 
wrap.  In  contrast  to  Fig.  3(a)  riddled  with  the  dark  fringes,  the 
retrieved  phase  map  Fig.  5(a)  does  not  present  any  dark  fringes 
at  all.  This  observation  is  also  verified  by  comparing  the  phase 
profile  Fig.  5(b)  to  the  intensity  profile  Fig.  3(b).  The  up-down 
bipolar  overshootings  appearing  in  Fig.  3(b)  disappear  in  the 
phase  profile  Fig.  5(b),  where  only  sharp  downward  peaks 
present  for  depicting  the  large  phase  shifts  at  the  bubble  rims. 
Therefore,  the  wrap’s  phase  map  Fig.  5(a)  clarifies  the  nature 
of  the  dark  fringes  in  Fig.  3(a)  such  that  these  dark  fringes  do 
not  indicate  any  stmctural  break  in  the  wrap  sample. 

IV.  DISCUSSION  AND  CONCLUSIONS 

This  study  demonstrates  that  the  robustness  against  noises 
of  a  phase  retrieval  method  is  critical  for  the  qualities  of  the 
retrieved  phase  map.  Specifically,  the  striking  differences 
between  the  retrieved  phase  maps  in  Figs.  4(a)  and  4(b), 
and  Fig.  5(a)  underline  the  importance  of  removing  the 
zero-frequency  singularity  that  is  intrinsic  to  the  TIE-based 
method.  While  the  TIE-based  phase  retrieval  method  is  compu¬ 
tational  effective  for  the  cases  with  strong  phase  contrast  effects 
and  low  noise  levels,"’'^  it  completely  failed  in  retrieving  the 
phase  map  of  the  bubble  wrap,  as  is  shown  in  Figs.  4(a)  and 
4(b).  This  is  because  the  noise  problem  is  especially  challeng¬ 
ing  for  the  case  at  hand,  as  the  bubble  wrap  has  very  low 
attenuation  contrast  and  only  moderate  x-ray  exposures  were 
applied  in  bubble  image  acquisitions.  From  mathematical  view¬ 
point,  the  TIE-based  phase  retrieval  formula  Eq.  (2)  is  ill-posed, 
since  it  involves  the  inverse  Laplacian  operator  that  owns 
a  zero-frequency  singularity.  Beyond  the  mathematics  formal¬ 
ity,  this  zero-frequency  singularity  reflects  the  fact  that  the 
phase  contrast  projection  is  insensitive  to  the  slow  variations  of 
x-ray  phase  shifts.  Seeking  to  recover  the  phase  shifts  by  com¬ 
paring  the  two  projection  images  contaminated  with  noise,  the 
TIE-based  method  has  amplified  the  noise  in  retrieving  the 
slowly  varying  phase  shifts  and  ruined  the  phase  retrieval  com¬ 
pletely.  The  Tikhonov  regularization  mends  the  singularity 
problem  by  replacing  the  inverse  Laplacian  operator  in 


Eq.  (2)  by  a  pseudodifferential  operator  V^/[(V^)^+a^]  with  a 
regularization  parameter  o:2>0.  As  is  shown  in  Fig.  4(b),  this 
regularization  made  some  but  only  little  improvement  com¬ 
pared  to  Fig.  4(a)  as  the  bubbles  start  to  appear  in  the  phase 
map  but  are  still  hardly  visible.  From  the  mathematical  formu¬ 
lation  of  the  Tikhonov  regularization  technique,  it  is  clear  that 
Tikhonov  regularization  just  limits  the  overall  magnitude  of  the 
errors  caused  by  the  noise  in  the  retrieved  phase  map,  but  the 
details  of  the  original  phase  map  could  still  get  lost  in  the  re¬ 
trieval.^^  The  AP-based  method  mends  the  singularity  problem 
by  utilizing  the  correlations  between  x-ray  attenuation  and  x- 
ray  phase  shifts  generated  by  tissues  or  materials,  as  is  shown 
in  Eq.  (8).  In  a  sense,  the  AP-based  method  utilizes  the  object’s 
x-ray  attenuation  to  get  rid  of  the  low-frequency  singularity  and 
the  associated  noise  amplification.  This  physics-motivated  reg¬ 
ularization  scheme  for  the  AP-based  method  is  expected  to  be 
more  effective  than  Tikhonov  regularization  in  phase  retrievals. 
The  greatly  improved  quality  of  the  retrieved  phase  map  shown 
in  Fig.  5(a),  as  compared  to  Figs.  4(a)  and  4(b),  validates  this 
expectation.  Hence,  the  performance  differences  of  these  phase 
retrieval  methods  underscore  the  importance  of  developing  an 
effective  method  to  remove  the  singularity  that  is  intrinsic  to 
the  TIE-based  method.  After  all,  the  high  noise  levels  in  the 
acquired  images  still  take  tolls  on  the  quality  of  the  phase  map 
Eig.  5(a)  retrieved  the  AP-method,  where  the  noise  are  quite 
visible.  While  increasing  radiation  doses  used  in  the  projections 
is  a  possible  solution  for  improving  phase  retrievals,  but  the 
radiation  dose  constraints  are  stringent  in  many  applications.  In 
clinical  applications,  it  is  imperative  to  limit  and  reduce  radia¬ 
tion  doses  involved  in  the  imaging.  Therefore,  more  research  is 
needed  for  developing  ever-improved  phase  retrieval  methods 
for  future  clinical  applications  of  phase  contrast  imaging. 

Other  factors,  which  are  less  critical  to  the  robustness  of 
phase  retrieval  but  may  affect  the  accuracies  of  retrieved 
phase  maps,  include  the  detector  calibrations  such  as  the  flat- 
field  and  gain  corrections,  and  the  ways  to  incorporate  the 
spectral  averaging  effects  of  polychromatic  x-rays.  In  fact, 
the  residual  background  nonuniformity  in  the  acquired 
images  [Figs.  2(a)  and  3(a)]  also  caused  variations  in  the 
background  phase  values  in  Fig.  5  and  reduced  their  accura¬ 
cies.  A  careful  detector  calibration  in  future  experiments 
should  avoid  this  kind  of  problem.  In  addition,  note  that  x- 
ray  wavelength  enters  as  an  important  parameter  in  Eq.  (2) 
for  the  TIE-based  method,  and  in  the  flow  chart  in  Fig.  I  for 
the  AP-based  method.  The  average  x-ray  wavelength  used  in 
the  phase  retrievals  should  represent  the  wavelength’s  linear 
and  nonliner  effects  averaged  over  the  exiting  x-rays  spec¬ 
trum  and  the  detector’s  spectral  response.  For  the  TIE-based 
method,  several  works  discussed  the  ways  of  performing  the 
spectral  averages  over  polychromatic  x-rays. We  did  not 
use  these  techniques  in  this  work,  because  the  TIE-based 
method  in  any  way  failed  the  wrap’s  phase  retrieval.  In  this 
work,  we  simply  use  the  estimated  average  photon  energy  of 
the  incident  x-ray  to  incorporating  the  spectral  averaging.  In 
the  future,  we  may  develop  more  elaborated  ways  for  incor¬ 
porating  the  spectral  averaging  effects.  Finally,  we  mention 
that  one  could  as  well  employ  a  special  phase  retrieval 
method  for  the  simple  samples  such  as  the  bubble  wrap.  This 


Medical  Physics,  Vol.  38,  No.  9,  September  2011 


5080 


Yan,  Wu,  and  Liu:  Robustness  of  phase  retrieval  methods 


5080 


special  method  works  specifically  for  the  single-material  ho¬ 
mogeneous  samples,  as  long  as  the  linear  attenuation  coeffi¬ 
cient  and  refractive  index  of  this  material  are  provided.  In 
this  special  method,  the  phase  map  of  a  single-material  sam¬ 
ple  can  be  retrieved  from  just  a  single  phase  contrast  image 
of  the  sample. The  air-bubble  wrap  can  be  approximately 
treated  as  a  single-material  sample  as  the  air  in  bubbles  con¬ 
tributes  negligibly  to  x-ray  attenuation  and  phase  shifts.  In 
this  work,  we  do  not  compare  this  special  method  for  single¬ 
material  samples  to  the  TIE-based  and  AP-based  methods, 
since  both  the  TIE-based  and  AP-based  methods  are  the  gen¬ 
eral  phase  retrieval  methods  applicable  for  any  general 
objects  of  inhomogeneous  materials. 

In  summary,  in  this  work,  we  compared  the  robustness  of 
two  phase  retrieval  methods,  the  TIE-based  method  and  the 
AP-based  method,  by  analyzing  the  retrieved  phase  maps 
from  the  experimental  projection  images  of  an  air-bubble 
wrap.  We  showed  that  the  TIE-based  method,  regardless  if 
the  Tikhonov  regularization  is  used,  failed  in  retrieving  the 
wrap’s  phase  map.  In  contrast,  in  the  wrap  phase  map 
retrieved  by  using  the  AP-based  method  bubbles  are  clearly 
recovered.  The  retrieved  phase  values  with  this  method  are 
reasonably  close  to  the  estimate  based  on  the  thickness- 
based  measurement.  The  stark  performance  differences  of 
the  two  methods  are  rooted  in  their  different  techniques 
employed  to  deal  with  the  singularity  problem.  This  compar¬ 
ison  shows  that  the  conventional  TIE-based  phase  retrieval 
method,  regardless  of  using  Tikhonov  regularization  or  not, 
is  unstable  against  the  noise  in  the  wrap’s  projection  images, 
while  the  AP-based  phase  retrieval  method  is  shown  in  these 
experiments  to  be  superior  to  the  TIE-based  method  for  the 
robustness  in  performing  the  phase  retrieval. 
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In  the  inline  phase  contrast  X-ray  tomography  the  reconstructed  apparent  linear  attenuation  coefficient 
values  may  be  greatly  larger  than  sample’s  linear  attenuation  coefficients  or  even  be  negative.  In  this 
work  we  present  a  general  formula  to  quantitatively  relate  the  apparent  linear  attenuation  coefficient 
values  in  cone-beam  phase  contrast  tomography  to  sample’s  linear  attenuation  coefficients  and  refractive 
indices.  This  formula  overcomes  the  gross  inaccuracy  of  the  existing  formula  in  the  literature  in  analyzing 
high-resolution  phase  contrast  tomography,  and  it  will  be  useful  for  correctly  interpreting  and  quantify¬ 
ing  the  apparent  linear  attenuation  coefficients  in  cone-beam  X-ray  phase  contrast  tomography. 
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For  many  applications  the  X-ray  attenuation  contrast,  such  as 
mass-density  differences  between  different  material  components. 
Is  too  small  to  generate  enough  image  contrast  with  conventional 
X-ray  tomography.  Taking  advantage  of  highly  sensitive  X-ray 
phase  shifts  generated  from  X-ray  coherent  scattering,  phase  con¬ 
trast  X-ray  tomography  finds  broad  applications  in  material  sci¬ 
ence  and  engineering,  biomedical  imaging  and  many  other 
scientific  fields  [1-3].  In  phase  contrast  X-ray  tomography  one  first 
makes  phase  contrast  angular  projections,  in  which  the  refraction 
and  diffraction  of  phase  shifted  X-rays  form  the  edge-enhancing 
interference  fringes  at  interfaces  of  different  material  components. 
From  these  angular  projections  3-D  tomograms  are  then  recon¬ 
structed  using  the  standard  tomography  reconstruction  methods 
such  as  the  filtered  backprojectlon  methods.  The  reconstructed 
phase-sensitive  tomograms  exhibit  enhanced  interfaces  and  layers 
of  different  materials  and  tissues  [1-3].  However,  these  tomograms 
are  not  the  maps  of  the  linear  attenuation  coefficients  (LACs)  of  the 
sample,  rather  they  are  the  maps  of  the  reconstructed  apparent 
LACs,  which  may  present  anomalously  large  or  even  negative 
apparent  LAC  values  at  Interfaces  between  different  materials 
and  tissues.  These  “artifacts”  may  cause  faulty  Interpretation  of 
sample  structures,  and  Impede  even  qualitative  characterizations 
for  tissues  and  materials  [l,2,4-6].  As  a  result,  for  example,  the 
apparent  LAC  values  of  small  bronchi  were  found  lower  than  that 
of  air  in  phase  contrast  tomography  experiments  [2].  In  order  to 
correctly  Interpret  and  quantify  the  mixed  contrast  exhibited  In 
phase  contrast  tomography,  it  is  necessary  to  find  out  quantitative 
relations  between  the  apparent  LAC,  on  the  one  hand,  and  the  ac¬ 
tual  LAC  and  refractive  index,  on  the  other.  The  previous  studies 


*  Corresponding  author.  Tel.:  +1  205  975  8135;  fax:  +1  205  975  4679. 
E-mail  addresses:  ayan@uabmc.edu  (A.  Yan),  xwu@uabmc.edu  (X.  Wu). 

0168-583X/$  -  see  front  matter  ©  2011  Elsevier  B.V.  All  rights  reserved. 
doi:l  0.1 016/j.nimb.201 1 .05.01 8 


found  that  the  apparent  LAC  in  X-ray  phase  contrast 

tomography  is  given  by  [7-9]; 

A‘„(ro)  =  /t(fo)-|.  g  +  ^  +  |ijA,e.(ro).  (1) 

In  this  equation  the  first  additive  term  ^((ro)denotes  the  sample’s 
actual  LAC  value  at  position  Pq,  and  the  second  additive  term  is  pro¬ 
portional  to  the  three-dimensional  Laplacian  of  the  sample’s  refrac¬ 
tive  Index  decrement  Arefr(ro)  [7-9[.  in  Eq.  (1)  R2  denotes  the 
sample-detector  distance  and  M  the  magnification  factor  employed 
in  angular  projections.  In  the  literature  this  second  term  is  used  to 
explain  the  edge  enhancement  and  the  anomalous  values  of  large 
apparent  LAC  [7-9[. 

However,  in  high-resolution  phase  contrast  tomography  with 
micrometer-scale  resolutions  the  applicability  of  Eq.  (1)  is  called 
into  question.  Researchers  recently  found  that  Eq.  (1)  is  invalid  in 
some  synchrotron  radiation  microtomography  experiments  [10]. 
In  fact,  in  the  derivation  of  Eq.  (1)  the  phase-sensitive  projections 
were  modeled  based  on  the  so-called  Transport  of  Intensity  Equa¬ 
tion  (TIE)  [7-9].  The  TIE  describes  how  the  phase-shifted  X-rays  re¬ 
fract  and  diffract  to  form  interference  fringes  at  boundaries  of 
material  components  [11,12].  However,  we  noted  that  TIE  can  only 
describe  the  moderate-resolution  phase-sensitive  projections  with 
TtAR2/(4Ma^)  «  1,  where  a  denotes  the  finest  size  resolved.  In 
cases  of  high-resolution  projections  with  micrometer-scale  resolu¬ 
tions,  the  fine  size  resolved  could  be  so  small  such  that 
nAR2/4Ma^  ~  1  or  larger,  then  TIE  becomes  Inapplicable  and  Eq. 
(1)  cannot  quantify  the  apparent  LAC  values  observed  in  high- 
resolution  phase  contrast  tomography. 

To  overcome  the  limitations  of  Eq.  (1),  we  derived  a  novel  and 
general  formula  of  apparent  LAC  values  In  phase  contrast  X-ray 
tomography.  Instead  of  using  the  TIE,  in  the  derivation  we 
employed  the  general  phase-sensitive  projection  formula  that  Is 
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applicable  for  hlgh-resolutlon  projections  as  well  [13-14],  We  then 
applied  the  standard  Feldkamp-Davis-Kresss  (FDK)  cone-beam 
reconstruction  algorithm  to  the  ray-sums  of  phase-sensitive  angu¬ 
lar  projections  [15],  We  found  that  the  apparent  LAC  value 
fiapparent{ro)  is  giveit  by 

,  f  aR2  ^  471 

f^appare„t('-o)  =  COS  ^ 

xsin(^V^o)A,efr(?o),  (2) 

where  the  three-dimensional  Laplaclan  operator  V3J,  =  {d^/dx^+ 
/dx\  +  /dx\).  The  operators  sin  (^{xk2/47tM)V3D)  and  cos({xR2 

/AnM)V\g)  in  above  equation  are  the  pseudo-differential  operators. 
The  action  of  a  pseudo-differential  operator  such  as  slnjcVjo)  on  a 
function  g(ro)  is  defined  as  sinjcV^ojgjfo)  =  /  /  exp  (i2ji 

(r-  to)  ■/)  ■  sin  (^-4n'^cp^g(f)d^fd^f.  Note  that  Eq.  (2)  is  indeed  a 

generalization  of  Eq.  (1).  For  experiments  with  moderate-resolu¬ 
tions  or  short  sample-detector  distances  such  that 
n2R2/4Ma^  «  1,  we  then  have  cos{{4R2/4nM)\/lu)  «  1  and 
sin((xR2/47tM)V3[,)  «  (/.R2/47cM)V3o.  Thereby  the  first  term  on 
the  RHS  of  Eq.  (2)  reduces  to  the  actual  LAC,  the  second  term  to  a 
scaled  Laplacian  of  Arefrjto)  and  Eq.  (2)  recovers  Eq.  (1),  as  is 
expected. 

In  order  to  validate  Eq.  (2),  we  performed  numerical  simula¬ 
tions.  In  the  simulations  the  X-ray  source  was  a  point  source  of 
60  keV  X-rays.  The  LACs  and  refractive  indices  of  the  numerical 
phantom  were  constructed  as  the  sums  of  the  LACs  and  refractive 
indices  of  six  spheres  of  light  elements  with  different  radii  and  cen¬ 
ter  positions.  For  each  sphere  the  mass  densities  were  designed  to 
decrease  radially  and  diminish  rapidly  at  the  sphere’s  boundary  for 
generating  phase  contrast.  Specifically  we  set  /i{ro)  =  0.215x 
ELi/n(?o)(l/cm)  and  Arefr(ro)  =  (7.548  x  10“®)  x  ELi/n(?o), 
where  /„{ro)  =  sjRl-(fo-  r„f/4R„  If  |{ro-r„)|^R„  and 
f„(ro)  =  0  otherwise.  Here  r„  and  R„  (n  =  1,  2 . 6)  denote  the  cen¬ 

ter  position  and  radius  of  the  respective  sphere.  In  this  way  the 
LACs  and  refractive  indices  change  continuously  without  any  steps 


Fig.  1.  Reconstructed  apparent  LAC  map  of  a  phantom  slice  from  the  phase- 
sensitive  cone-beam  tomography  with  mixed  contrast. 


but  vary  rapidly  at  the  sphere  boundaries.  Hence  strong  phase  con¬ 
trast  would  manifest  around  the  boundaries  of  the  spheres.  The 
detector  has  512  x  512  pixels  with  a  sampling  pitch  of  2  pm.  In 
the  simulated  scan  the  source  traced  a  circle  orbit,  and  the  source 
and  detector  rotated  together  over  an  angular  range  of  360°  in  0.5° 
steps.  In  order  to  simulate  X-ray  Fresnel  diffractions  Involved  in 
the  angular  projections,  one  needs  to  compute  X-ray  attenuation 
values  and  phase  shifts  of  the  exit  rays  in  each  of  the  projections. 
In  general  one  needs  to  use  the  multlsllce  propagation  method  to 
account  for  X-ray  Fresnel  diffractions  inside  the  sample  [16],  How¬ 
ever,  the  settings  described  above  satisfy  the  so-called  thin  object 
condition  T  <  p^/2xo,  where  T  denotes  the  sample’s  thickness,  p  Is 
the  resolution  and  xq  the  wavelength  of  60  keV  X-ray  [16],  This 


Fig.  2.  Profiles  of  the  apparent  LACs  along  the  vertical  line  passing  the  center  of  the  central  sphere,  (a)  Profile  measured  from  Fig.  1.  (b)  Profile  calculated  based  on  Eq.  (1). 
(c)  Profile  calculated  based  on  the  general  formula  Eq.  (2). 
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Table  1 

A  comparison  of  the  fringe  peak-peak  differences  in  Fig.  2(a)-(c). 


Fringe  set 

Fringe  set 

Fringe  set 

Fringe  set 

Fringe  set 

Fringe  set 

Fringe  set 

Fringe  set 

1 

2 

3 

4 

5 

6 

7 

8 

Peak-peak  difference  of  apparent  LACs  in  Fig.  2(a)  (1/cm) 

4.648 

4.979 

6.098 

8.860 

8.860 

6.098 

4.979 

4.648 

Error  of  calculated  peak-peak  difference  based  on  Eq.  (1) 

429 

346 

199 

136 

136 

199 

346 

429 

(Fig.  2(b))  (%) 

Error  of  calculated  peak-peak  difference  based  on  Eq.  (2) 

15.0 

17.5 

13.0 

6.9 

6.9 

13.0 

17.5 

15.0 

(Fig.  2(c))  (%) 


condition  allows  us  to  compute  the  X-ray  attenuation  values  and 
phase  shifts  of  the  exit  X-rays  as  the  ray  integrals  of  the  sample’s 
LACs  and  refractive  index  decrements.  In  this  simulation  a  ray¬ 
tracing  algorithm  was  used  to  compute  the  ray  integrals  such  as 
f  fidl  and  /  Arefrdi  in  each  of  the  angular  projections.  The  exit 
X-ray  attenuation  and  phase  shift  were  determined  as  = 
exp  (-/ /idl)  and  cp  =  -(Inj).)  f  Arefrdi,  respectively  [11-14].  Using 
the  exit  X-ray  attenuation  and  phase  shift  from  a  given  projection, 
we  then  simulated  the  free-space  Fresnel  diffractions  of  the  exit  X- 
rays  propagating  from  the  sample  to  the  detector  with  a  setting  of 
R2  =  1  m  and  M  =  2  [16].  The  Fresnel  diffraction  intensities  in  all  the 
angular  projections  were  log-transformed  to  calculate  the  ray- 
sums  needed  for  tomographic  image  reconstruction.  We  applied 
the  standard  FDK  algorithm  for  image  reconstruction  [15].  Fig.  1 
shows  a  reconstructed  slice  of  apparent  LAC  values  in  gray  scale. 
This  slice  shows  the  edge-enhancing  dark-bright  fringes  caused 
by  phase  contrast  at  the  interfaces  between  the  spheres.  Fig.  2(a) 
shows  the  profile  of  the  apparent  LAC  values  along  the  vertical  cen¬ 
tral  line  of  Fig.  1.  The  ordinate  denotes  apparent  LAC  value  in  (1/ 
cm)  along  the  vertical  line,  while  the  abscissa  denotes  pixel  loca¬ 
tion  along  the  line.  For  the  sake  of  comparison  Fig.  2(b)  and 
Fig.  2(c)  show  the  profiles  of  the  apparent  LAC  values  along  the  ver¬ 
tical  line,  which  were  computed  by  using  Eq.  (1 )  and  Eq.  (2)  respec¬ 
tively.  Note  that  the  ordinate  scale  in  Fig.  2(b)  is  three  times  larger 
than  that  of  Fig.  2(a)  and  (c).  To  facilitate  the  comparisons  of  the 
three  profiles,  we  measured  the  fringe’s  peak-to-peak  differences 
(the  changes  between  peak  and  trough)  at  eight  fringe  sets  for  each 
of  the  three  profiles  shown  in  Fig.  2(a)-(c).  Table  1  lists  the  mea¬ 
sured  fringe  peak-to-peak  differences  of  apparent  LACs  in 
Fig.  2(a),  which  is  obtained  from  Fig.  1,  the  phantom’s  phase  con¬ 
trast  tomogram  simulated  by  using  the  Fresnel  diffraction  and 
FDK-reconstruction.  Compared  to  these  eight  peak-to-peak  differ¬ 
ences  in  Fig.  2(a),  the  table  lists  as  well  the  errors  of  the  calculated 
peak-peak  differences  based  on  Eq.  (1)  (Eig.  2(b))  and  Eq.  (2) 
(Eig.  2(c))  respectively.  Erom  these  comparisons,  it  is  clear  that 
the  apparent  LAC  values  calculated  from  Eq.  (1 )  are  grossly  inaccu¬ 
rate  with  the  errors  amounting  to  several  hundred  percent,  while 
the  apparent  LAC  values  calculated  by  means  of  Eq.  (2)  achieve  rea¬ 
sonable  quantitative  accuracies  as  is  shown  in  Table  1.  In  addition, 
the  profile  Eig.  2(c)  shows  two  small  bumps  at  the  center  of  the 
profile.  We  found  that  these  variations  of  apparent  LACs  are  arti¬ 
facts  caused  by  omitting  the  small  nonlinear  terms  in  phase  shifts 
in  the  derivation  of  Eq.  (2).  The  omission  of  these  nonlinear  terms 
mainly  affects  the  accuracies  in  computing  the  side-lobes  of  phase 
contrast  fringes,  thereby  causes  these  two  small  bumps,  which  are 
associated  with  the  distorted  fringe  side-lobes  near  the  right 
boundary  of  the  left  off-center  sphere  in  the  phantom  (Eig.  1). 

In  summaiy,  we  present  a  general  formula  (Eq.  (2))  to  quantita¬ 
tively  relate  the  apparent  LAC  values  in  cone-beam  phase  contrast 


tomography  to  sample’s  LAC  and  refractive  index  values.  This 
formula  overcomes  the  gross  inaccuracy  of  the  existing  formula 
in  the  literature  in  analyzing  high-resolution  phase  contrast 
tomography.  We  believe  that  this  general  formula  will  be  useful 
for  designing  phase  contrast  tomography  experiments,  and  for  cor¬ 
rectly  interpreting  and  quantifying  the  apparent  LAC  values  in 
cone-beam  X-ray  phase  contrast  tomography.  For  example,  Eq. 
(2)  can  be  used  for  computing  the  edge-enhancement  profiles  for 
analyzing  presampling  modulation  transfer  functions  in  synchro¬ 
tron  radiation  microtomography  [10]. 
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